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ABSTRACT 

We have recently completed a search of 6 years of archival BATSE data for gamma- 
ray bursts (GRBs) that were too faint to activate the real-time burst detection system 
running onboard the spacecraft. These "non-triggered" bursts can be combined with 
the "triggered" bursts detected onboard to produce a GRB intensity distribution that 
reaches peak fluxes a factor of ~ 2 lower than could be studied previously. The value 
of the (V/V max ) statistic (in Euclidean space) for the bursts we detect is 0.177 ± 0.006. 
This surprisingly low value is obtained because we detected very few bursts on the 4.096 
s and 8.192 s time scales (where most bursts have their highest signal-to-noise ratio) that 
were not already detected on the 1.024 s time scale. If allowance is made for a power- 
law distribution of intrinsic peak luminosities, the extended peak flux distribution is 
consistent with models in which the redshift distribution of the gamma-ray burst rate 
approximately traces the star formation history of the Universe. We argue that this 
class of models is preferred over those in which the burst rate is independent of redshift. 
We use the peak flux distribution to derive a limit of 10% (99% confidence) on the 
fraction of the total burst rate that could be contributed by a spatially homogeneous (in 
Euclidean space) subpopulation of burst sources, such as type Ib/c supernovae. These 
results lend support to the conclusions of previous studies predicting that relatively few 
faint "classical" GRBs will be found below the BATSE onboard detection threshold. 
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Introduction 

The origin of some, and possibly all, gamma-ray 
buiists (GRBs) at cosmological distances has been 



firmly established with the identification of X-ray op- 



tical, and radio after 


glows (Costa et al. 1997 




Van 


Paradijs et al. 1997 




Frail et al. 1997 


) and the sub- 



sequent measurement of cosmological redshifts for 
at least four o f the optical afterglows and/or their 



1998a 



host galaxies (Mctzger et al. 1997; Kulkarni et al 



Djorgovski et al. 1999j ; pjorgovski et al. 1998| ) 



Tht: objects responsible tor producing the majority ot 
GRBs, the gamma-ray bursters themselves, have yet 
to be understood, however. To obtain an understand- 
ing of the spatial distribution of sources and the dis- 
tribution of their burst luminosities is a crucial step 
towards identifying the physical processes that pro- 
duce GRBs. 

Before the rapid follow-up of GRB afterglows was 
made possible by the BeppoSAX satellite, the only 
way to test hypotheses about the spatial and lumi- 
nosity distributions was to fit parametric models to 
the measured characteristics of the bursts themselves. 
For this purpose the distribution of GRB intensities 



was used (see, for example, 


?enimore et al. 1993; 


Rut- 


.edge, Hui, & Lewin 1995; 


Fcnimore & Bloom 1995|; 


Cohen & Piran 1995; 


Hakkila et al. 1996 


; and refer- 



ences therein). The effects of cosmological time dila- 
tion on the time profiles of bright versus faint bursts 
were also studied ( Norris et al. 1995| ) . Since optical 
spectroscopic redshifts are so far associated with only 
four (possibly five) bursts^, number counts as a func- 
tion of intensity remain an important tool for explor- 
ing the possible spatial and luminosity distributions 
of GRBs. 



Several recent papers (Totani 1997, 1998; 


Wijers 


3t al. 1998 




Krumholz, Thorsctt, & Harrison 199£|; 


Mao & Mo, 1998) have used the observed GRB in- 



tensity distributions to investigate the possibility that 
the redshift distribution of gamma-ray bursters traces 
the global star formation history of the Universe. The 
motivation for this hypothesis is a collection of theo- 
retical models in which GRBs are produced by stel- 
lar objects that evolve from their formation to their 
bursting phase on a time scale of ~ 100 Myr or less. 
This group of models includes the merging of a neu- 



lr The proposed association of GRB 980425 with SN 1998bw (z 
= 0.008; Galama et al. 1998) may indicate a separate class of 
GRBs (Bloom et al. 1998). We will therefore consider that 
event separately (see section 3.2). 



tron star with another neutron star or a black hole, 
the collapse of a massive star, and the collapse of a 
Chandrasekhar-mass white dwarf (see Wijers et al 



1998 for references). In these scenarios, the cosmo- 
logical redshift distribution of the GRB rate should 
approximately follow the redshift distribution of the 
formation rate of stellar objects; in other words, the 
GRB rate should trace the global star formation his- 
tory of the Universe. This hypothesis appears to solve 
some puzzling aspec ts of the observations, such as the 
•'no h ost" problem ( Schacfcr et al. 1997| ; Wijers et al 



The star formation rate (SFR) as a function of 
redshift has been studied by Lilly et al. (1996), Fall, 
Chariot, & Pei (1996), Madau, Pozzetti, & Dickinson 
(1998), and Hughes et al. (1998). The principal re- 
sult of these studies is that the SFR was substantially 
higher in the past. Between the present and z ss 1 the 
SFR increases by a factor of ~ 10; it peaks somewhere 
in the range z ~ 1 to z rs 3; and it decreases to a rate 
comparable to the present by z « 4-5 (this last point 
remains uncertain). 

Totani (1997), Wijers et al. (1998), Krumholz et al. 
(1998), and Mao & Mo (1998) all find that the GRB 
peak flux number counts can accommodate the hy- 
pothesis that the GRB rate follows the SFR. Among 
the important conclusions that these authors derive 
from this interpretation of the data are the follow- 
ing: 1) that the faintest gamma-ray bursts observed 
with the Burst and Transient Source Experiment 
(BATSE) onboard the Compton Gamma Ray Obser- 
vatory (CGRO) are already produced at redshifts of 
z w 3 to z 6 (Wijers et al. 1998; but see Section 0); 



and 2) that more sensitive experiments are unlikely to 
discover large numbers of faint GRBs (of the "classi- 
cal" kind that are detected with current instruments) 
below the BATSE onboard detection threshold. The 
latter conclusion has important implications for the 
design and operation of future GRB detectors, which 
will test the behavior of GRB number counts at in- 
tensities well below the BATSE threshold. 

We have recently completed a search of 6 years 
of archival data from BATSE for GRBs and other 
transients that did not activate the real-time burst 
detection system (or "trigger") running onboard the 
spacecraft. A GRB or other transient may fail to ac- 
tivate the BATSE onboard burst trigger for any of 
several reasons. The burst may be too faint to exceed 
the onboard detection threshold, it may occur while 
the onboard trigger is disabled for technical reasons, 



2 



it may occur while the onboard trigger is optimized 
for detecting non-GRB phenomena, or it may arti- 



fici dly raise the onboard background estimate and 



be mistaken for a below-threshold event. Our search 
of the archival data is sensitive to GRBs with peak 
fluxes (measured over 1.024 s in the 50-300 keV en- 
ergy range) that are a factor of ~ 2 lower than can be 
detected with the onboard trigger in its nominal con- 
figuration. Thus our search constitutes an experiment 
that is ~ 2 times more sensitive than those reported 



in the BATSE catalogs ( 


Fishman et al. 1994; Meegan 


et il. 1996; 


Paciesas et al. 1999 




Meegan et al. 1998 





In this paper we present results regarding the peak 
flux distribution of the GRBs detected with our "off- 
line" search of archival data. In section || we sum- 
marize some important aspects of our off-line search 
and we discuss the (V/V max ) statistic for the bursts 
we detected. We show that surprisingly few bursts 
are found on the 4.096 s and 8.192 s time scales that 
were not already detected on the 1.024 s time scale. 
In section |^ we fit parametric cosmological models 
to the observed differential peak flux distribution to 
compare scenarios in which the GRB rate follows the 
SFR with the model in which the co-moving GRB 
rate is independent of redshift. We also examine the 
possibility that a homogeneous (in Euclidean space) 
population of bursting objects could be contributing 
to the observed sample of GRBs. In section^ we show 
how our results provide two independent arguments 
that favor models in which the GRB rate follows the 
SFR over models in which the GRB rate is indepen- 
dent of redshift. 

2. The Search for Non- Triggered GRBs 



The details of our off-line search of the BATSE 
data are discussed in Kommers et al. (1997). We have 
merely extended the search from covering 345 days 
of the mission to covering 2200 days. We have also 
made minor modifications to our peak flux estimation 
procedure in order to secure better relative calibration 
between our peak fluxes and those in the 4B catalog 
( Paciesas et al. 1999| ) . The extended catalog of non- 
triggered events will be provided and discussed in the 
Non- Triggered Supplement to the BATSE Gamma- 
Ray Burst Catalogs (Kommers et al. 1999b). Here 
we address only those aspects of the search that are 
rel evant to the GRB intensity distribution analysis. — 



with 1.024 s time resolution (the data type desig- 
nated "DISCLA" in the flight software; Fishman et 
1989). These data are searched for statistically 



al 



significant count rate increases to identify candidate 
burst events. The many candidate events ("off-line 
triggers") are then visually inspected to separate as- 
tronomically interesting transients from instrumental 
and terrestrial effects. To be considered a GRB, a 
candidate must exhibit significant signal in the 50-300 
keV range (DISCLA channels 2 and 3) and it must 
lack any characteristics that would associate it with 
a solar flare, Earth magnetospheric particle precipi- 
tation, or other non-GRB origin. Since the DISCLA 
data are (nearly) continuously recorded, our search 
detects some bursts that already activated the on- 
board burst trigger; we call these events "onboard- 
triggered bursts." Bursts that were detected exclu- 
sively by our search of archival data are called "non- 
triggered bursts." 

In addition to searching at the 1.024 s time reso- 
lution of the DISCLA files, we also search the data 
binned at 4.096 s and 8.192 s time resolution. The 
longer time bins provide greater sensitivity to faint 
bursts that have durations longer than ~ 4 s or ~ 8 
s. The specific time profile of each burst determines 
which of these three time scales is the most sensitive. 
For this reason the searches on each time scale should 
be considered separate experiments. 

Our search covers 1.33 x 10 8 s of archival data span- 
ning the time from 1991 December 9 to 1997 Decem- 
ber 16. In these data we detected 2265 GRBs, of 
which 1392 activated the onboard burst trigger and 
873 did not. We will refer to these 2265 GRBs as the 
"off-line GRB sample." During the same time period, 
the onboard burst trigger detected 1815 GRBs. The 
1815 — 1392 = 423 bursts that were detected by the 
onboard burst trigger but that were not detected by 
our search either occurred during gaps in the archival 
DISCLA data or had durations much less than the 
1.024 s time resolution (so they did not achieve ade- 
quate statistical significance in the archival data). 

Note that because the best time resolution available 
to our retrospective search is 1.024 s ; a ^ results in 
this paper pertain to bursts with durations longer than 
about 1 s. Thus, the population of "short" (duration 
less than ~ 2 s) bursts that contributes to the bi- 



modal GRB duration distribution ( Kouvcliotou et al 



— use — Lire data fioiu — the Large Area Delec- 

tors that provide count rates in 4 energy channels 



1993| ) is not well represented in the off-line sample. 



An estimate for the fraction of bursts that our search 
is likely to miss because of our time resolution can be 
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obtained from the 4B catalog. Although 21% of GRBs 
for which both durations and flucnces were available 
had T 90 < 1.024 s, only 7% had both T 90 < 1.024 s 
and fluences too small to create adequate statistical 



significance in the 1.024 s data (Paciesas et al. 1999). 

For each of the 873 non-triggered GRBs we have es- 
timated a peak flux in the 50-300 keV range based on 
the time bin with the most counts above background. 
For 1288 of the 1392 onboard-triggered GRBs, we 
used the peak fluxes from the current BATSE GRB 
catalog (Mecgan et al. 1998). For the remaining 104 



onboard-triggered bursts, peak fluxes were not avail- 
able from the current burst catalog; we estimated 
peak fluxes for them using our own techniques as we 
did for the non-triggered bursts. 

Since the onboard trigger criteria were changed for 
a variety of reasons during the time spanned by our 
search, we adopt for the nominal onboard detection 
threshold the value 0.3 ph cm" 2 s" 1 in the 50-300 keV 
range. At this peak flux the on board trigger efficiency 
is w 0.5 ( Paciesas et al. 1999 ). With this estimate, 
551 of our 873 non-triggered bursts were below the 
nominal onboard detection threshold. The rest were 
not detected onboard for the reasons cited previously. 

2.1. Trigger Efficiency 

To determine the peak flux threshold of the off-line 
GRB sample, the trigger efficiency E\ (P) of our off- 
line search has been calculated using the techniques 
described in Kommers et al. (1997). This quantity is 
the probability that a burst that occupies exactly one 
1.024 s time bin with a peak flux P will be detected 
by the off-line search algorithm. E\(P) is well rep- 
resented (within the uncertainties of the calculation 
owing to variations in the background rates) by the 
following function: 



•erf (-3.125 +16.677P)], (1) 



where erf(x) is the standard error function and P is 
given in units of photons cm~ 2 s -1 in the 50-300 keV 
band. This equation is plotted as the dashed line in 
Figure [ij Error bars on the grid points of the calcula- 
tion (diamonds) represent the sample standard devia- 
tion of the calculated probabilities owing to variations 
in the background rates. For comparison, the BATSE 
trigger efficiency from the 4B catalog ( Paciesas et al 



be detected, however. This is because many GRBs in 
our sample last longer than 1.024 s; therefore, these 
bursts have more than one statistical chance to be 
included in the sample. 

Suppose the peak of a burst occupies N time bins, 
so that the burst has effectively N statistical chances 
to be detected. Then the probability that the burst 
is detected can be approximated as unity minus the 
probability that the burst fails to be detected in all 
N trials: 



E N {P) =1- [l-Ei(P)} 



N 



(2) 



Since the number of chances N is not known for GRBs 
a priori, the actual probability of detection E(P) is 
obtained by marginalizing En(P) over the distribu- 
tion of N for bursts with peak fluxes near P: 



E(P) = 



J2h(N,P)E N (P) 

ZHn,p) ■ 



(3) 



Our estimate for h(N, P), the histogram of the various 
integer values of N for bursts with peak fluxes near P, 
was obtained from the detected sample of bursts by 
counting, for each burst, the number of time bins with 
count rates that were within one standard deviation 
of the peak count rate. For purposes of illustration, 
Figure [| shows the histogram of N for bursts with 
peak fluxes in the range 0.1-0.4 ph cm~ 2 s _1 . The 
resulting function E{P) is well represented (to within 
the uncertainties of the calculation) by the formula 

E(P) = - [1 + erf (-4.801 + 29.868P)] . (4) 

This equation expresses our best estimate of the trig- 
ger efficiency of our off-line search on the 1.024 s time 
scale. It is plotted as the solid line in Figure [l]. The 
efficiency of our search falls below 0.5 at a peak flux 
of 0.16 ph cm~ 2 s _1 . 

If we had not made some correction for the effect 
of time profiles on the single-time-bin burst detec- 
tion probabilities, we would have substantially un- 
derestimated our trigger efficiency near the detection 
threshold (~ 0.2 ph cm~ 2 s _1 ). We note that this 
type of correction to the single-time-bin trigger effi- 
ciency should also be applied when using the trigger 
efficiencies given in the IB, 2B, 3B, and 4B catalogs 
(Fishman et al. 1994; Meegan et al. 1996; Paciesas 



19£ 9 ) is plotted as the dotted line (grid points are et al. 

indicated by open squares). Equation |l| tends to un- in't Zand & 

derestimate the probability that a typical GRB will man (1995) 



1999). Similar considerations are addressed by 
Fenimore (1994) and Loredo & Wasser- 
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2.2. (C mm /C max ) 3/2 Distribution 




0.0 L 



0.1 0.2 0.3 0.4 

Peak Flux (ph crrT 2 s _1 50-300 keV) 

Fig. 1. — Trigger efficiency for our off-line search. 
The grid points of the calculations are plotted as 
individual symbols. Error bars represent the stan- 
dard deviations of the calculated probabilities owing 
to variations in the background rates. The dashed 
line (equation 1) shows the probability that a burst 
occupying a single time bin is detected by our search. 
The solid line (equation 4) shows the marginal prob- 
ability that a burst is detected by our search, given 
that some bursts longer than 1.024 s have more than 
one statistical chance to be detected. For comparison, 
the dotted line shows the trigger efficiency from the 
4B catalog; no uncertainties are available for the grid 
points (squares). 



As successively more sensitive instruments have 
been used to produce GRB catalogs, it has been cus- 
tomary to give the val ue of the (V/V m ax) statistic 
for th e detected bursts ( Schmidt, Higdon, fc Huetei 
198S ) . For photon counting experiments like BATSE, 
it is not strictly (V/V max ) that is typically calcu- 
lated, but rather ((C m in/C m ax) 3/2 }, where C min is 
the threshold count rate and C m ax is the maximum 
count rate measured during the burst. The depar- 
ture of ((C m in/C max ) 3 / 2 ) from the value of h ex- 
pected for a population of bursters distributed homo- 
geneously in Euclidean space (with a well-behaved, 
but otherwise arbitrary luminosity distribution) has 
been firmly established ( Mccgan et al. 1992| ; Mee- 
gan et al. 1996Q. Since the discovery that most 



GRBs originate at cosmological distances, the quan- 
tity {{C minj Cmax) 3 1 2 ) can no longer be interpreted 
as {V/V m ax)- Nevertheless, it is useful to compare 
the values of ({C m in/Cmax) 3 ^ 2 } obtained by succes- 
sively more sensitive experiments, including the value 
obtained for the bursts detected with our search. 

Table | lists various missions and the values they 
obtained for ((C m in/Cmax) 3 ^ 2 ) ■ The trend towards 
lower values of ({C m in/ C m ax) 3 ^ 2 ) with more sensitive 
experiments indicates that increasing the accessible 
survey volume by decreasing the flux threshold does 
not lead to the detection of large numbers of faint 
bursts. 

The value of ((C* rm „/C ma:E ) 3/2 ) for the 2265 GRBs 
detected by our search|] is 0.177 ± 0.006. This is 
the lowest value ever obtained for a sample of GRBs. 
The cumulative distribution of (C m inl C m ax) 3 ^ 2 for 
our GRBs is shown in Figure [|. The flattening of 
this curve in the range 0.5 < [C m inl C m ax) 3 ^ 2 < 1-0 
shows that over 90% of the GRBs we detect are above 
threshold (on at least one of the 3 time scales) by a 
factor of at least (0.5)~ 3 / 2 = 1.6. 

The reason for this low value of ((C m in /C ma x) 3 ^ 2 ) 
is the fact that most of the bursts we detected had 
their maximum signal-to-noise ratios on the 4.096 s 
and 8.192 s time scales, yet surprisingly few bursts 
were detected only on these longer time scales. For 
each burst we compute the values of (C m i n /Cmax) 3 ^ 2 
on each of the 3 time scales. The largest of the three 
values for each burst is used in taking the average. In 



2 This value supersedes the ones given in Kommers et al. (1996, 
1997, 1998), which are incorrect due to a programming error. 
An erratum has been submitted (Kommers et al. 1999a). 
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Fig. 2. — Histogram of N, the number of time bins 
within one standard deviation of the peak count rate, 
for bursts with peak fluxes in the range 0.1-0.4 ph 




Fig. 3. — Cumulative distribution of (C m in/C m ax) 3 ^ 2 
for the off-line GRB sample. The dramatic flatten- 
ing of the curve above (C m i n /C max ) 3 ^ 2 — 0.5 shows 
that few of the GRBs detected in our search are just 
barely above the detection threshold on all 3 time 
scales (1.024 s, 4.096 s, and 8.192 s). 



Euclidean space this corresponds to taking for each 
burst the smallest value of (V/V max ). Sinc e 72.0% 
of the bursts we detected have Tqq durations (Koshut 



et al. 1996) longer than 8 s, we expect the average 
((Cmin/Crnax) 3 ^ 2 ) to be dominated by values mea- 
sured on the 8.192 s time scale. 

In fact, the average ((C mm /C' max ) 3/2 ) = 0.177 ± 
0.006 includes 520 values measured on the 1.024 s 
time scale, 491 values measured on the 4.096 s time 
scale, and 1254 values measured on the 8.192 s time 
scale. Yet only 105 bursts were detected exclusively 
on either of the 4.096 or 8.192 s time scales (or both). 
Many of the bursts that are barely above the detec- 
tion threshold on the 1.024 s time scale are well above 
the detection threshold on the longer time scales. 
Thus very few bursts are found to be just barely 
above our detection threshold on all 3 time scales, and 
this accounts for the low value of ({C m i n /C max ) 3 ^ 2 ) . 
Restricting our calculation to use only count rates 
measured on the 1.024 s time scale (and bursts de- 
tected on the 1.024 s time scale) gives a larger value, 
{{C min /C max ) 3 / 2 ) = 0.247 ± 0.006. 

Roughly, the 4.096 s search should be ~ 2 times 
more sensitive than the 1.024 s search for bursts that 
maintain their peak flux for at least ~ 4 s, and the 
8.192 s search should be yet more sensitive. Therefore 
our lack of GRB detections exclusively on the longer 
time scales indicates either 1) a substantial paucity of 
faint, long bursts below the threshold of our 1.024 s 
search, or 2) that during our visual inspection of the 
off-line triggers we have tended to classify a substan- 
tial number of faint, long GRBs as other (non-GRB) 
phenomena. We feel that both alternatives must be 
present at some level. 

A review of the non-GRB off-line triggers suggests 
that events resembling faint, long GRBs that illumi- 
nate the same detectors as a known, bright, variable 
X-ray source are more likely to be attributed to vari- 
ability from the X-ray source than to be classified as 
GRBs. There is also a tendency to classify bursts 
that have directions consistent with the Sun as solar 
flares. A secondary evaluation of the event classifi- 
cations suggests that between 50 and 200 (this range 
represents the central 90% confidence interval) , with a 
most likely value of 86, GRBs have been misclassified 
in this way. The corresponding "loss rate" is between 
2% and 8% (most likely 4%) of the total 2265 bursts 
in the off-line sample. This is not enough to fully ex- 
plain, as experimental error, the dearth of faint, long 
bursts below our 1.024 s threshold. 
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2.3. Peak Fluxes 

Detailed comparisons of cosmological models with 
the data require intensity distributions in physical 
units. We have chosen to do the analysis in terms 
of the burst rate as a function of peak photon flux 
measured over 1.024 s in the energy range 50-300 keV. 
Compared with the fiuence (total energy per unit area 
deposited in the detector by the burst) we prefer peak 
photon flux for the purposes of intensity analysis. The 
peak photon flux can be obtained more reliably from 
the raw count data and it is more directly related to 
our ability to detect bursts. 

Of the 2265 GRBs detected by our search, we chose 
to include in our peak flux analysis only those that 
were detected on the 1.024 s time scale, so that equa- 
tion ^ gives the detection efficiency. We also chose to 
use only those bursts with peak fluxes in the range 
0.18-20.0 ph cm -2 s _1 . The lower limit ensures that 
the off-line trigger efficiency exceeds 0.8 for the range 
of intensities used in the analysis, and the upper limit 
excludes very bright bursts which are too rare to pro- 
vide adequate counting statistics in narrow peak flux 
bins. With these cuts on the data, we are left with 
1998 peak flux measurements. To fit the differential 
intensity distribution, we bin the 1998 bursts into 25 
peak flux intervals that were chosen to be approx- 
imately evenly spaced in the logarithm of P. The 
spacing is A log P w 0.05 in the range 0.18 < P < 1.0, 
AlogP w 0.1 in the range 1.0 < P < 7.9, and there 
is a final broad bin for the range 7.9 < P < 20.0. 
Uncertainties in the number of bursts AN b s in each 
bin are taken to be ±\/ AN Q b s . The burst rate is com- 
puted by dividing the number of bursts in each bin 
by the live time of the search (1.33 x 10 s s = 4.21 yr) 
and the mean solid angle visible to the BATSE detec- 
tors (0.67 x 4-71"). Table || gives the peak flux intervals, 
number of bursts, and burst rate for each bin. 

3. Cosmological Model Comparison 

Many investigators, in scores of papers, have shown 
the consistency of the GRB peak flux distribution 
with various cosmological models (see, for exam- 



pie, 


Wijers et al. 1998|; Lorcdo & Wasserman 1998]; 


Hakkilaet al. 1996; 


Horacket al. 1996 


; Rutlcdge, Hui, 


k Lewin 1995; 


Fcnimore & Bloom 1995|; and refer- 



ences therein). As shown in the previous section, the 
off-line GRB sample extends the observed GRB in- 
tensity distribution to peak fluxes that are lower by a 
factor of ~ 2 than could be studied previously. While 



it is unlikely that a factor of ~ 2 will yield stringent 
new model constraints, it remains of interest to note a 
few cosmological models that provide good fits to the 
extended GRB peak flux distribution. These can be 
used to set limits on the rate of GRBs that may come 
from a nearby, spatially homogeneous subpopulation 
of burst sources. 

3.1. Purely Cosmological Models 

To limit the number of free parameters that must 
be considered, our choice of cosmological world model 
is the Einstein-de Sitter model (O = 1, A = 0, (fa = h\ 
Weinberg 1972). This cosmology has been used by 



many other investigators so it allows easy comparison 
of results. Where needed, we assume a Hubble con- 
stant H = 70 /170 km s _1 Mpc -1 . We also assume 
that bursters are distributed isotropically, so the only 
interesting parameter in the burster spatial (redshift) 
distribution is the radial coordinate r(z) from Earth. 
The following derivation of the expected peak flux 
distributions follows the discussions in Fenimore & 
Bloom (1995) and Loredo & Wasserman (1997). 

In general the rate of bursts R per unit interval in 
peak flux P observable at Earth is given by 



dR 
dP 



— = dL dz 



8 2 R 
dL dz 



6(P-$(L,z)), (5) 



where L is the equivalent isotropic peak luminosity of 
the burst at the source, z is the redshift parameter, 
d 2 R/dL dz is the rate of bursts per unit L per unit 
redshift interval, 5(x) is the Dirac delta function, and 
Q(L,z) is the peak photon flux measured at Earth 
for a burst with peak luminosity L located at redshift 
z. We will assume that the redshift and luminosity 
distributions are independent, so that the burst rate 
as a function of L and z is given by 



d 2 R 
dL dz 



4:7TCR 



H 



Hz) 



(l + z) 2 



(G) 



where i?o is an overall normalization, tp(L) is the dis- 
tribution of burst luminosities (normalized to unity), 
p{z) is the distribution of the co-moving burst rate 
as a function of redshift (normalized to unity on the 
interval < z < 10), and r(z) = (2c/H ){\ + z - 
\/l + z)/(l + z) is the co-moving radial coordinate. 

The peak flux $(£, z) observed at Earth in the 50- 
300 keV energy range, where the BATSE burst trigger 
is sensitive, depends on the intrinsic spectrum of the 
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GRB. We write it as 



*(L,z) = 



L K(z) 



47r (1 + z) r 2 (z) 



(7) 



The spectral correction function K(z) depends on the 
shape of the burst photon energy spectrum at the 
source. The observed GRBs have a variety of spec- 
tral shapes, and in the cosmological scenario these 
observed spectra have been redshifted according to 
the (unknown) redshifts of the sources. 

To account for the spectral variety of GRBs we use 
the spectral fits of Band et al. (1993). To account for 
the unknown redshift factors for these spectra, we use 
the procedure described in Fenimorc & Bloom (1995). 
The peak fluxes of the bursts for which Band et al. 
(1993) derived spectral fits are used in conjunction 
with the cosmological model under consideration to 
self-consistently estimate the redshift factors for the 
fitted spectra. We assume that the ith burst fitted 
by Band et al. (1993) has exactly the mean intrin- 
sic peak luminosity in the cosmological model being 
considered: Li = JdL' L'tp(L'), where the shape of 
ip(L) depends on the parameters of the model lumi- 
nosity function. We then solve for the redshift Zi 
which the fitted burst i must have had to produce 
the peak flux listed for it in the current BATSE GRB 
catalog (Meegan et al. 1998). Fifty-one of the bursts 



fitted by Band et al. (1993) had peak fluxes available. 
For each of their spectral shapes 4>i{E) the spectral 
correction function takes the form 



Ki(z) = 



f 300(l+ Z )/(l+ Zi ) ,„ , ,„s 

J 50( i +z)/( i +Zi) ah <k[E) 
(1+*) J^l^dEEME) 



(8) 



The integrals in the denominator and numerator con- 
vert the model parameter L, which represents the 
peak luminosity in the 30-2000 keV range at the 
source, to the observed photons cm -2 s _1 in the 50- 
300 keV band at Earth. The burst rate expected in 
the BATSE band pass for the ith spectral shape <j>i(E) 
is then (from equation 0) 



(dR\ 
fdz 



16tt 2 cR 



i) 



P(z) r\z) 



0) 



(1+z) VT+^ K 



,/. ( ^ + r 2 M P \ 



The limits on the integral are determined by the range 
of z for which -0(47r(l + z)r 2 (z)P/Ki(z)) is non-zero 
at the given P. 

To estimate the observed distribution of bursts, 
which includes a variety of spectral shapes, we average 



Equation || over the 51 spectral correction functions 
Ki(z). This procedure is equivalent to marginaliz- 
ing the unknown spectral parameters of the observed 
bursts (i.e., those in the off-line sample) to obtain the 
posterior rate distribution. The 51 spectra from Band 
et al. (1993) are furnishing estimates of the prior dis- 
tributions of the spectral parameters. The expecta- 
tion value of the observed burst rate for peak fluxes 
between Pi and Pi is then 



AR(P U P 2 ) = 
j£dPE(P) (§), 



(10) 



where (dR/dP) is the mean rate estimated from the 
51 observed spectra and E(P) is the detection effi- 
ciency. 

The use of the Band et al. spectra increases the 
computational cost of the rate model by a factor 
of ~ 50 over using a single "universal" burst spec- 
trum. We found that a simple power-law form for the 
GRB photon energy spectrum — as has been used by 
many previous studies — predicts significantly differ- 
ent burst rates at low peak fluxes than does equation 
|Io| . Since we are interested in the behavior of the 
burst rate at low peak fluxes, we felt that the anal- 
ysis based on the full 51 Band et al. spectra would 
be more reliable. Similar conclusions are reached by 
Fenimore & Bloom (1995) and Mallozzi, Pendleton, 
k Paciesas (1996). 

For comparison with the results of previous stud- 
ies, we chose two forms for the luminosity distribu- 
tion. The first is a monoluminous (standard candle) 
distribution. The second is a truncated power-law, 



I Llog(L max /L„ 
1p{L) = i (1-/3) L-0 







= 1 



(11) 



with ip(L) = if L < L m i n or L > L max . The nor- 
malization factors ensure that JdL tp(L) = 1. 

The standard candle distribution, though useful for 
comparison with other results, is ruled out by the ob- 
served peak fluxes of the four bursts for which as- 
sociated optical redshifts have been measured. For 
GRBs 970508, 971214, 980613, and 980703 the in- 
ferred equivalent isotropic peak luminosities in the 
30-2000 keV energy range are given in Table |. To 
calculate each of these peak luminosities, we have 
used the observed 50-300 keV peak flux (on the 1.024 
s time scale) in combination with the observed red- 
shift to find the expectation value of the intrinsic lu- 
minosity averaged over the 51 Band et al. spectra. 
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This procedure is the one used in our modeling, so it 



consistent with recent results from SCUBA (Hughes 



was used on these four bursts also, to facilitate com- 
parisons with the models (see section |4|) . The peak lu- 
minosities estimated here are somewhat higher by fac- 
tors of ~ 3 to ~ 6 than those reported elsewhere (e.g., 
Krumholz et al. 1998). This is because the spectral 
shapes fitted by Band et al. (1993) generally become 
steeper at high energies, so a source at high redshift 
must be more luminous to produce the flux observed 
at Earth than it would have to be if the spectrum did 
not fall off so rapidly at higher energies. These dif- 
ferences illustrate the importance of using the most 
realistic spectral models available rather than simple 
power-laws when analyzing the GRB intensity distri- 
bution. 

A variety of spatial, or rather redshift, distribu- 
tions for the bursters have been used in previous stud- 
ies of the GRB intensity distribution. With up to 4 
free parameters already incorporated into our burst 
rate models (the overall normalization i?o, and the 
parameters of the power-law luminosity function j3, 
L m in, and L max ) there is little hope of constraining 
any additional free parameters in the redshift dis- 
tribution. Here we explore 3 specific models of the 
redshift distribution that contain no free parameters. 
The two physical scenarios we examine are 1) that 
the co-moving burst rate is independent of redshift 
between z = and z — 10, and 2) that the co-moving 
GRB rate is proportional to the star formation rate 
(SFR). 

For the GRB rate model that is independent of 
redshift, p(z) = 0.1 for < z < 10 and p{z) = 
for z > 10. We refer to this redshift distribution as 
"model Dl." 

For the case where the burst rate follows the star 
formation history of the Universe, we use two slightly 
different parameterizations of the SFR. The first is the 
SFR deduced from the rest-frame ultraviolet luminos- 
ity density, with the functional form given in footnote 
1 of Madau, Delia Valle, & Panagia (1998). In this 
estimation the SFR peaks around z =1—1.5. A SFR 
of roughly this form has been used by several previ- 



et al. 1998) which are not susceptible to the same 



ous studies of the GRB intensity distribution (Totani 
19S0; |Wijers et al. 1998fc Krumholz et al. 1998). We 
refer to this redshift distribution as "model D2." 

The second SFR parameterization assumes that 
the SFR — and thus the GRB rate — tracks the to- 
tal output of radio-lou d AGN. In this scenario the 
SFR peaks at z =2-3 ([Hughes et al. 1998j |Dunlop 
19S(8|). This form of the SFR appears to be more 



problems of dust obscuration as the determination by 
Madau et al. (1998). The specific functional form 
we use is a best-fit analytic model to points mea- 
sured by hand from Figure 6 of Hughes et al. (1998): 
p(z) cx 0.00360 + 0.0108 exp(2.76z - 0.573z 2 ). This 
approximation appears to be accurate to within 5% 
for the redshift range 1 < z < 4. (At lower and higher 
redshifts the formula likely underestimates the actual 
rate of star formation; but this is no great concern 
as it is the redshift of the peak SFR that is of pri- 
mary interest.) We refer to this redshift distribution 
as "model D3." 

With choices for ^){L) and p(z) as discussed, we 
fit equation [Io| to the data in Table || by minimizing 
the x 2 statistic. In all cases, we found that the pa- 
rameter L max was not well constrained: variations in 
L m ax did not change the minimum \ 2 by a significant 
amount. The (mathematical) reason for this is that 
the integrand in equation ^ is a decreasing function 
of z for plausible values of 0, so that varying the up- 
per limit {zmax corresponding to L max for the given 
P) causes only small changes in the value of the in- 
tegral. Accordingly, all the results reported here set 



lOOOZv 



The free parameters are thus Rq 



and Lq in the cases of the standard candle models, 
and i?o , /3, and L m i n in the cases of the power-law lu- 
minosity distribution models. The results of the fits 
are listed in Table [| and Table ^. Uncertainties on 
the fitted parameters correspond to 68% confidence 
limits for 2 (A^ 2 = 2.3) or 3 (A% 2 = 3.5) interesting 
parameters, respectively ( Avni 1976| ). 

Model Dl (constant burst rate density as a func- 
tion of redshift) produces an acceptable fit for the 
standard candle luminosity distribution. The proba- 
bility of getting x 2 > 32.3 for 23 degrees of freedom 
is 0.094. Adding one more free parameter (/?) for 
the power-law luminosity distribution produces an in- 
significant change in the minimum x 2 ■ Furthermore, 
the high value of f3 in the best-fit power-law distribu- 
tion indicates a very narrow range of peak luminosi- 
ties. 

Model D2 (burst rate density follows the SFR as 
determined by Madau et al. 1998) produces a formally 
unacceptable fit in the monoluminous case. But it 
achieves an excellent fit (x 2 per degree of freedom = 
0.81) for the power-law luminosity distribution. The 
F-test estimates a probability of 1.5 x 10~ 7 that the 
improvement in x 2 is due to chance, justifying the 
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inclusion of the additional parameter in the power- 
law luminosity function model. The value of in this 
model is remarkably well-constrained. If the other fit 
parameters are regarded as "uninteresting" then the 
90% (Ax 2 = 2.7) an d 99% (A X 2 = 6.6) confidence 
intervals flAvni 1976| ) on (3 are 2.0-2.3 and 1.8-2.6, 
respectively. 

Model D3 (burst rate density follows the output 
of radio-loud AGN) produces formally acceptable fits 
with both the standard candle and power-law lumi- 
nosity distributions. The power-law luminosity dis- 
tribution achieves a significantly lower y 2 , however. 
The -F-test estimate of the probability that the im- 
provement is due to chance is 1.5 x 10~ 3 . 

Figure |] plots the differential peak flux distribu- 
tions for the best-fit models with power-law luminos- 
ity distributions. For all three best-fit models the 
value of {(Pmin/P) 3 ^ 2 ) is consistent with the value of 
({C m in I 'CmaxY' 2 ) measured for the sample (see sec- 
tion ^j). Extrapolating the best-fit models to peak 
fluxes lower than those included in our data shows 
very different behaviors. Model Dl (dot-dashed lines) 
predicts a dramatically higher burst rate at low peak 
fluxes than do models D2 (solid line) and D3 (dashed 
line). 

In each model the best-fit parameters for the power- 
law luminosity function yield our best estimate of the 
parameters of the intrinsic distribution of GRB peak 
luminosities. The distribution of peak luminosities of 
the observed bursts is different, however, because the 
most luminous bursts are sampled from a much larger 
volume than than are the least luminous bursts. Even 
though high luminosity bursts are infrequent, the ge- 
ometrical advantage of sampling them from a larger 
volume means that they will be over-represented in 
a sample of bursts observed over a fixed time inter- 
val. The distribution of peak luminosities for the ob- 
served bursts is the "effective luminosity function" 
(see Loredo & Wasserman 1997 for further discus- 



sion). For the best-fit parameters of model Dl the 
effective luminosity function is a power-law that is 
less steep than that of the intrinsic luminosity func- 
tion. We find (3^g = 2.8 for the effective luminosity 
function versus (3 = 4.6 for the intrinsic one. The 
power-law slopes of the effective luminosity functions 
for models D2 and D3 are = 1.6 and = 1.9, 
respectively. 

Similarly, the distribution of the GRB rate as a 
function of redshift for the observed bursts is not iden- 
tical to the intrinsic redshift distribution given by p(z) 




constant 

SFR {Madau et al. 1998) 
SFR (Hughes et al- 1998) 



-0.5 0.0 
00,10 F (ph cm 



0.5 1,0 

' 50-300 keV) 



Fig. 4. — Best-fit cosmological models with power- 
law luminosity distributions. Units of R are bursts 
yr _1 sr _1 , and those of P are ph cm~ 2 s _1 in 50-300 
keV. The dot-dashed line corresponds to model Dl 
(co-moving burst rate is independent of redshift). The 
solid line shows model D2 (burst rate follows the rest- 
frame ultraviolet luminosity density) and the dashed 
line shows model D3 (burst rate follows the output of 
radio-loud AGN). Measured rates are shown with la 
vertical error bars; horizontal error bars indicate the 
bin widths. The best-fit model curves displayed here 
have not been corrected for detection efficiency. 
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(see Loredo & Wasserman 1997). Figure || shows the 
effective redshift distributions for the best-fit mod- 
els Dl, D2, and D3. In all 3 models the effective 
redshift distribution cuts off at a lower redshift than 
does the corresponding intrinsic redshift distribution. 
The mean redshifts of the observed bursts in the best- 



fit models Dl, D2, and D3 are 



\D1 



0.86, (z. 



D2 



1.4, and (z) D3 — 1.9, respectively. The maximum 
redshifts of the bursts in models Dl, D2, and D3 are 
not precisely determined, but the observed redshift 
distributions cut off around z^* x ~ 1.5, z™ x ~ 2, 



and 



r D3 



3, respectively (see Figure 



It is interesting to compare the rate of bursts that 
are seen with the off-line search to the total rate of 
bursts that occur in the Universe, subject to the cos- 
mological rate models we are considering. The frac- 
tion of bursts that are detected with BATSE is given 
by the integral of the effective redshift distribution 
for the off-line search (over the range < z < 10) 
divided by the integral of the effective redshift distri- 
bution that would be visible to a "perfect detector" 
which can detect infinitely faint bursts, multiplied by 
the live-time (0.70) and sky exposure (0.67) fractions. 
In the best-fit models Dl, D2, and D3, it follows that 
the off-line search detects ~ 14%, ~ 30%, and ~ 37%, 
respectively, of the bursts that occur in the Universe. 
Of course, these estimates can apply only to bursts 
with energy spectra and durations of the kind that 
are accessible to the off-line search. If BATSE had 
100% live-time and no Earth blockage, then in model 
Dl 30% of the bursts that occur could be detected 
with the off-line search. In models D2 and D3, 65% 
and 78%, respectively, of the bursts that occur could 
be detected with the off-line search. 

It is customary to quote the rate of GRBs as the 
co-moving rate per unit volume at z = 0, a quantity 



often denoted by po with units of Gpc 3 yr (Feni- 



mo f-e fc Bloom 1995| , |Wijers et al. 1998[ ). The model 
parameter i?o is related to po by po — 47ri?o/o(0). Ta- 
ble H lists the values of po corresponding to the best- 
fit values of i?o for the standard-candle GRB models. 
Table does the same for the GRB models with a 
power-law luminosity function. This burst rate can 
be converted into an event rate per "typical" galaxy 
using the space density of such galaxies. Here, a "typ- 
ical" galaxy is taken to be one with a luminosity L* 
equal to the characteristic luminosity of galaxies in 
the Schechter function ( Schechter 1976 ). Loveday et 
al. (1992) report the mean space density of L* galax- 
ies to be (4.8 ± 0.6) 10~ 3 h? m Mpc~ 3 . With this 



m 
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Fig. 5. — Redshift distribution of the burst rate for 
observed bursts (i.e., the effective differential burst 
rate as function of redshift). The distributions are 
normalized so that the integral of the distribution vis- 
ible to a "perfect" detector with no sensitivity limit is 
unity. The dot-dashed line corresponds to a constant 
burst rate as a function of redshift (model Dl). The 
solid line corresponds to the model (D2) where the 
burst rate traces the SFR as determined by Madau 
et al. (1998). The dashed line corresponds to the 
model (D3) where the burst rate traces the output 
of radio-loud AGN (Hughes et al. 1998). In model Dl 
the off-line search detects a smaller percentage of the 
bursts that occur than it does in models D2 and D3. 
The redshifts associated with GRBs 970508, 971214, 
980613, and 980703 are marked with asterisks. 
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conversion, the co-moving rate of GRBs at z = can 
be expressed in GEM (Galactic Events per Million 



years), which is the rate of GRBs in an L* galaxy. 
This quantity is also listed in Tables |6| and |7|. 

3.2. Limits on a Possible Homogeneous Sub- 
population 

The discovery of the unusual Type Ib/c supernova 
SN 1998bw in the X-ray error box of GRB 980425 



duced the GRB ( 


Galama et al. 1998; 


Iwamoto et al. 


1998 




Kulkarni ct al. 1998b|). It has been suggested 



that such events (the supernova- GRBs, o r "S-GRBs") 



may constitute a subclass of all GRBs ( Bloom et al 
1998). In this subsection we discuss what fraction 



of all GRBs could belong to such a subclass assuming 
that the remaining bursts come from the "reasonable" 
cosmological scenarios discussed previously. 

If the inferred peak luminosity of GRB 980425 
(assuming a distance corresponding to the redshift 
z = 0.0085 of SN 1998bw; |Tinney et al. 1998|) is 



typical of S-GRBs, then the bursts in this subclass 
are detectable within a volume of radius ~ 100 Mpc 
( Bloom ct al. 199"§| ). Within this volume, we assume 
the spatial distribution of Type Ib/c SNe to be ap- 
proximately homogeneous. Thus the cumulative in- 
tensity distribution of S-GRBs can be expected to 
follow a —3/2 power-law. This conclusion follows for 
any well-behaved distribution of intrinsic luminosities 
as long the spatial distribution does not deviate from 
homogeneity within the volume sampled by our de- 
tectors. Since the observed intensity distribution of 
all GRBs deviates strongly from the —3/2 power-law, 
it can be used to set an upper limit on the fraction of 
all GRBs that can come from a subclass which obeys 
the —3/2 power-law. In this respect, the faint end of 
the peak flux distribution (as explored by our off-line 
search) provides the most stringent constraints. 

A model-independent limit on the fraction of bursts 
that might come from a nearby homogeneous popula- 
tion can be obtained from the histogram of (C m i n / C max ) 3 ^ 2 
for the observed bursts. Any homogeneous subpopu- 
lation is expected to contribute a constant number of 
bursts to each bin in the histogram. Thus the bin with 
the fewest number of bursts sets an upper limit on 
the total number of bursts that could come from the 
subpopulation. The bursts detected with our search 
ha ve already been used to set an upper limit of 5- 



an upper limit on the fraction of bursts that could 
be associated with nearby Type Ib/c SNe (Kippen et 
al. 1998). Though model- independent, this limit de- 



pends on how coarsely the histogram is binned. It also 
assumes than an arbitrary distribution of intensities 
is acceptable for bursts that do not come from the ho- 
mogeneous subpopulation. This is too much freedom, 
because physically plausible distributions occupy only 
a subset of all arbitrary intensity distributions. 

Here we assume that the bulk of GRBs come f rom 
the cosmological distributions discussed in section |3.1| . 
An upper limit on the rate of all GRBs that might 
come from Type Ib/c SNe (or any nearby homoge- 
neous distribution) is then fixed by determining the 
maximum rate of bursts that can come from a (differ- 
ential) distribution proportional to P~ 5 ' 2 before the 
model becomes inconsistent with the data. We thus 
fit a model of the form 



AR(P 1 ,P 2 ) = R H J^dP E(P) P- 6 ' 2 
J^dPE(P) ($). 



+ 



(12) 



The fractional burst rates corresponding to the 90% 
and 99% confidence upper limits on the normalization 
Rh in each model are given in Table |[ The upper 
limits were determined by finding the value of Rh 
for which A\ 2 = 2.7 and A% 2 = 6.6, respectively, 
when x 2 is minimized with respect to the other fit 
parameters ( Avni 1976 ). In all cases, only a modest 
fraction, 5-10%, of the observed GRBs could come 
from a homogeneous subpopulation (and thus from 
nearby SNe). These upper limits are comparable to 
the model-independent result found in the previous 
paragraph. They are slightly less constraining be- 
cause of the fact that our peak flux distribution refers 
only to the 1.024 s time scale, so it takes no account 
of the paucity of faint bursts found on the 4.096 and 



6% on the fraction of bursts that could come from 



8.192 s time scales (see section 2.2). 

These results were to be expected from the facts 
that 1) models Dl, D2, and D3 with power-law lu- 
minosity distributions already gave excellent fits to 
the data without the presence of the homogeneous 
(p -5 / 2 ) term, which is sharply peaked at low peak 
fluxes, and 2) the fractional uncertainties on the rates 
in each bin are on the order of 10%. The upper lim- 
its discussed here would be further reduced if a given 
GRB must exhibit certain characteristics (e.g., single- 
peaked time profile, lack of emission above 30 keV) 
in order to be considered a candidate S-GRB ( Bloom| 
et al. 19981) . Norris, Bonnell, & Watanabe (1998) 



a homogeneous subpopulation; this limit is therefore 



have found that only 0.25-0.5% of BATSE GRBs have 
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temporal and spectral characteristics similar to GRB 
980425. 

4. Discussion 

The GRB peak flux distribution alone (on the 
1.024 s time scale) only weakly distinguishes between 
the non-evolving model p{z) = constant (Dl) and the 
evolving models where p{z) is proportional to an esti- 
mate of the star formation history (D2 and D3). (In 
this section we restrict our attention to models that 
include a power-law distribution of intrinsic peak lu- 
minosities.) A similar conclusion has been reached 
previously by Krumholz et al. (1998), who analyze the 
BATSE catalog data and find that to reliably distin- 
guish the non-evolving and evolving models requires 
data from more sensitive GRB detectors and/or the 
measurement of more individual GRB redshifts. 

The off-line sample of GRBs constitutes a more 
sensitive experiment than the one analyzed by Krumholz 
et al. (1998). Here we argue that models similar (or 
identical) to D2 and D3, in which the GRB rate has 
a significant peak in the redshift range 1 < z < 3, 
are modestly preferred over the constant rate density 
model. Two independent lines of reasoning serve to 
denigrate the p(z) = constant ("non-evolving") mod- 
els in favor of the evolving ones. 

First, our search on the 1.024 s time scale can reach 
peak fluxes as low as 0.16 ph cm -2 s _1 in the 50- 
300 keV band (50% detection efficiency). But our 
searches on the 4.096 and 8.192 s time scales are sen- 
sitive to peak fluxes (averaged over the matching time 
scale) that are lower by factors of ~ 2 and ~ 2y/2, re- 
spectively, than the 1.024 s threshold. In fact, most 
of the bursts we detect have their highest signal-to- 
noise ratio in the 8.192 s search. Yet surprisingly few 
bursts are detected exclusively on the longer 4.096 s 
and 8.192 s time scales. This suggests that there are 
relatively few faint GRBs waiting to be detected by 
a search that is more sensitive than the one we car- 
ried out. In this respect, the evolving models (D2 and 
D3) appear to be more accurate. They predict that 
the number of bursts per logarithmic peak flux inter- 
val will level off towards lower peak fluxes, and may 
even start to decline (see Figure ||). The non-evolving 
model, on the other hand, predicts that the number 
of bursts observed per logarithmic peak flux interval 
will continue to increase towards lower peak fluxes. 

Since we have not derived peak fluxes on the 4.096 
and 8.192 s time scales in order to repeat the anal- 



ysis of section 3.1, we offer the following quanti- 
tative evidence that the paucity of bursts detected 
on the longer time scales favors the evolving mod- 
els D2 and D3 over the non-evolving model Dl. 
On the 1.024 s time scale, the measured value of 
{{C mm /C max ) 3 ' 2 ) = 0.247 ± 0.006 is trivially con- 
sistent with the values of {{Pmin/ P) 3 ^ 2 ) found for 
the best-fit models in section |3.1| . But the value 
of {{C mm /C max ) 3 / 2 ) = 0.177 ± 0.006 found for all 
bursts detected by our search contains information 
on the paucity of faint bursts on the 4.096 and 8.192 
s time scales. We can compare it with the value of 
{{Pmin/ P) 3 ^ 2 ) obtained by extrapolating the best-fit 
models of section 1A to the peak flux threshold asso- 
ciated with the 8.192 s search. Taking 0.18/(2^2) = 
0.06 ph cm -2 s _1 as the approximate P m in for the 
8.192 s search, we obtain the following values for 
{{P min /P) 3/2 )- 0.221 in model Dl, 0.169 in model 
D2, and 0.147 in model D3. Thus model D2 produces 
a value of {{P m in/ P) 3 ^ 2 ) that is the most consistent 
with the value of {{C min / C max ) 3 / 2 ) found for our full 
sample, and model Dl produces the most inconsistent 
value. 

Second, the inferred equivalent isotropic peak lu- 
minosities (in the 30-2000 keV range) of the 3 bursts 
for which associated redshifts have been measured 
can be compared with the effective luminosity dis- 
tributions of the best-fit models. The best-fit non- 
evolving model (Dl) predicts that 90% of all GRBs 
should come from the narrow range of intrinsic peak 
luminosities (0.29-0.66) lO 51 /^ 2 erg s _1 (a factor of 
~ 2). The range from which 90% of the observed 
GRBs in this model are drawn is somewhat broader, 
however: (0.29-1.5) 10 51 /if 2 erg s" 1 (a factor of ~ 
5). In contrast, the intrinsic luminosities inferred 
from the 3 bursts with associated redshift information 
span a much broader peak luminosity range, (0.6- 
37) 10 51 /if 2 erg s" 1 (a factor of ~ 62). Bursts with 
peak luminosities as high as those inferred for GRB 
971214 ((37 ± 16) 10 51 /if 2 erg s" 1 ) and GRB 980703 
((2.2±0.4) 10 51 h^Q erg s _1 ) are extremely rare events 
if the (non-evolving) model Dl is correct. On the 
other hand, the best-fit (evolving) models D2 and D3 
allow a much broader (and uniformly higher) range 
of luminosities. The best-fit model D2 predicts that 
90% of all GRBs are drawn from the intrinsic peak 
luminosity range (0.50-7.3) 10 51 /if 2 erg s _1 (a fac- 
tor of - 15) and that 90% of the observed GRBs 
are drawn from the range (0.53-72) 10 51 h^Q erg s _1 
(a factor of ~ 130). Likewise, the best-fit model D3 
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predicts that 90% of all GRBs are drawn from the 
range (1.5-9.3) 10 51 /if 2 erg s _1 (a factor of ~ 6) 
and 90% of the observed GRBs are drawn from the 
range (1.5-40) lO 51 /^ 2 erg s" 1 (a factor of - 25). In 
the context of models D2 and D3, the GRBs 970508, 
971214, 980613, and 980703 constitute a much more 
likely sample of detected bursts than in model Dl. 

The effective redshift distributions (the rates of ob- 
served bursts as a function of redshift) furnish an- 
other point of comparison with GRBs 970508, 971214, 
980613, 980703. As shown in Figure | the model Dl 
predicts a vanishingly small rate of observed bursts 
from the redshift z = 3. 418 measured for GRB 981214 
(Kulkarni ct al. 1998a). If model Dl were correct, 



then it would be remarkable that BATSE and Bep- 
poSAX detected such a rare burst: only 1 x 10 -6 
of the rate distribution comes from higher rcdshifts. 
On the other hand, the evolving models D2 and D3 
predict much higher rates of observed bursts from 
z = 3.418. Even in these models, however, such a 
high redshift is exceptional: in model D2 only 0.4% 
of the burst rate distribution lies beyond z = 3.4, 
and in model D3 only 1.8% does. Still, as in the case 
of the inferred luminosity distributions, the 3 bursts 
with associated redshifts are a much more likely sam- 
ple in models D2 and D3, where the GRB rate follows 
an estimate of the star formation history. 

These results certainly do not prove that the GRB 
rate traces the star formation rate. The peak flux dis- 
tributions alone fail to exclude the non-evolving rate 
density model (Dl) with high confidence, especially in 
view of the unknown cosmological parameters which 
can be varied to improve the fit. Furthermore, un- 
til many more redshifts are associated with specific 
GRBs and/or more sensitive GRB detectors go on- 
line, the data will not be able to distinguish qualita- 
tively similar SFR evolution models such as D2 and 
D3. Any evolution that specifies a significant peak 
in the burst rate in the redshift range 1 < z < 3 is 
likely to be consistent with current data. The rea- 
soning regarding the paucity of faint bursts detected 
only on the 4 s and 8 s time scales should be addressed 
more quantitatively in the context of the models; but 
this is difficult owing to our poor understanding of 
the diverse time profiles of GRBs and of the correla- 
tions between time profiles and peak fluxes. Finally, 
we have considered only 3 very specific "straw-man" 
models, and it may be that none of them are partic- 
ularly accurate representations of the true GRB rate 
density and peak luminosity distributions. For ex- 



ample, a previous episode of star formation at high 
redshift could contribute a hitherto undetected pop- 
ulation of very faint GRBs. 

Nevertheless, the results of our search appear to 
support the conclusions reached by Totani (1997), 
Wijers et al. (1998), and Krumholz et al. (1998), who 
showed that if the GRB rate traces the SFR then 
relatively few faint "classical" GRBs are to be found 
below the BATSE onboard detection threshold. This 
information should be useful to the designers and op- 
erators of future GRB detectors. 

In Figure ^ we plot the cumulative rate distribu- 
tion of the off-line sample of GRBs along with the 
best fit models and their extrapolations to lower and 
higher peak fluxes than are available in the data. The 
best fit models apparently predict too many bright 
bursts above the maximum peak flux that was used 
to fit the differential burst rate (20.0 ph cm -2 s _1 ). 
This discrepancy is not surprising because the fits did 
not include data at high peak fluxes. A Kolmogorov- 
Smirnov (K-S) test ( Press et al. 1995 ) indicates that 
the discrepancy is not significant if the data and best 
fit models are compared for all peak fluxes above 0.18 
ph cm~ 2 s _1 . (The probability of getting a larger 
value of the K-S statistic is 0.4 for model D2.) If the 
K-S test is restricted to the data above 10.0 ph cm~ 2 
s _1 , however, the discrepancy is significant (the prob- 
ability of getting a larger value of the K-S statictic is 
6 x 10~ 4 for model D2). This means that the extrap- 
olations of the best fit models are inaccurate at high 
peak fluxes. 

Future research should explore the peak flux dis- 
tribution over the widest possible range. For exam- 
ple, the GRB detector onboard Pioneer Venus Or- 
biter (PVO) operated for a much longer mission than 
has BATSE (so far) and it has more completely sam- 
pled the rate of very bright bursts. When the peak 
fluxes of the GRBs detected with PVO are calibrated 
to match the BATSE peak fluxes, it would be of in- 
terest to see if the parameters of the best-fit models 
found in this paper remain consistent with the num- 
ber counts when the very bright bursts are included. 
Figure || shows that the slopes of the best-fit mod- 
els approach the —3/2 slope rep orted to be consisten t 
with the brightest PVO bursts ( Fenimorc et al. 1993| ). 

The results of the best-fit models are otherwise 
generally consistent with previous studies of the BATSE 
data. In particular we find that while we cannot con- 
strain the full width of the power-law luminosity func- 
tion in any of the scenarios we considered, the best-fit 



14 



CL 
A 



cn 
o 




-2-10 1 2 

log 10 P (ph cm" 2 s _1 50-300 keV) 

Fig. 6. — Cumulative peak flux distributions for best- 
fit models. The units of R are bursts yr _1 sr _1 . The 
observed peak flux distribution for the off-line sample 
is shown as the solid histogram. The best-fit model 
Dl is shown as the dot-dashed line. Model D2 is 
shown as the solid line. Model D3 is shown as the 
dashed line. The discrepancy at high peak fluxes is 
not surprising because the fitting used only data from 
the interval 0.18-20.0 ph cm~ 2 s~ x , but the curves in 
the figure are intended to reflect true cumulative dis- 
tributions (with no lower or upper limit on peak flux). 
The best-fit models appear to predict too many bright 
bursts, though the discrepancy between the data and 
the models is not statistically significant unless the 
comparison is restricted to bright bursts only (data 
above 10.0 ph cm -2 s _1 ). This suggests that the off- 
line sample should be combined with data from longer 
missions, such as PVO and/or Ulysses to check the 
behavior of the burst rate at high peak fluxes. 



models yield intrinsic peak luminosity functions that 
contain 90% of all GRBs within a factor of 10-20. 
This result (or a similar one) has been previously ob- 
tained by Ulmer, Wijers, & Fenimore (1995), Woods 
& Loeb (1995), Hakkila et al. (1996), and Horack et 
al. (1996). We find that the peak luminosity distri- 
bution of the observed bursts, however, is wider: in 
model D2, 90% of the observed bursts come from a 
peak luminosity range that spans a factor of ~ 100 
or more. A similar result is discussed by Loredo & 
Wasserman (1998). 

The effective redshift distributions that we obtain 
(see Figure [5|) are reasonably consistent with those 
obtained by Krumholz et al. (1998) and Mao & Mo 
(1998) using the BATSE catalog data. They are also 
consistent with limits on the redshifts of GRB sources 
set by the non-detection of any gravitationally lcnsc d 
GRBs in the BATSE catalogs flMarani et al. 1998|) . 
Holz, Miller, & Quashnock (1998) derive upper lim- 
its of (z) = 2.3 (68% confidence) and (z) = 5.3 (95% 
confidence) for the average redshift of GRB sources 
in the Einstein-de Sitter cosmology. The values of (z) 
we obtain for our best fit models (see Section [O]) are 
all well within these limits. The effective redshift dis- 
tributions are also consistent with the disparity found 
by Norris et al. (1995) between the duration distribu- 
tions of bright and dim GRBs, which they interpret 
as the signature of cosmological time dilation for a 
GRB source distribution with the dimmest bursts at 
z « 2. 

On the other hand, our effective redshift distribu- 
tions are somewhat at odds with the conclusions of 
Wijers et al. (1998) that the faintest bursts observed 
with BATSE are at redshifts of 3 < z < 6. Part 
of the discrepancy is explained by the fact that rate 
model used by Wijers et al. (1998) does not correct 
for the cosmological time dilation of the co-moving 
burst rate. According to our results with models Dl, 
D2, and D3, no more than 0.002%, 1%, or 5% (re- 
spectively) of the bursts observed with our off-line 
sample are from redshifts larger than z = 3.0; and 
fewer than 8 x 10~ 7 %, 0.07%, and 0.3%, respectively, 
are at redshifts greater than z = 4.0. So, it is possible 
that GRBs may be produced at redshifts as high as 
z ss 6; but even model D3, which allows the highest 
redshifts, permits only a 7% chance that one or more 
such bursts have been observed among the 2265 in 
the off-line sample. 

To summarize, our off-line search of archival BATSE 
data has explored the distribution of GRB intensities 
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at peak fluxes below the onboard detection thresh- 
old. We find a paucity of faint bursts detected on the 
4.096 and 8.192 s time scales that were not already 
detected on the 1.024 s time scale. The differential in- 
tensity distribution is consistent with models in which 
the GRB rate traces the global star formation his- 
tory of the Universe; and it is marginally consistent 
with the model in which the GRB rate is independent 
of redshift. We argue that the models in which the 
GRB rate traces the star formation rate are never- 
theless preferred, based on the paucity of faint bursts 
detected exclusively on the 4.096 s and 8.192 s time 
scales, and on the comparison of the inferred effective 
luminosity and redshift distributions with the bursts 
for which redshifts have been measured. As an appli- 
cation of the off-line GRB intensity distribution, we 
set a limit of 10% (99% confidence) on the fractional 
rate of all GRBs that could belong to a homogeneous 
(in Euclidean space) subpopulation of burst sources 
(such as Type Ib/c supernovae). 
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Table 1 

Values of ((C mi „/C maa; ) 3/2 ) obtained by various GRB detectors. 



Detector/Mission 



Reference 



PVO 

Konus/ Venera 11 & 12 
KA/HEAO 1 
GBS/SMM 
GBD/ Ginga 
BATSE/ CGRO (3B) 
BATSE/ CGRO (off-line, 

1.024 s search only) 
BATSE/ CGRO (all off-line) 



0.46 ±0.02 
0.45 ±0.03 
0.40 ±0.08 
0.40 ±0.025 
0.35 ±0.035 
0.33 ±0.01 

0.247 ±0.006 
0.177 ±0.006 



Hartmann et al. 1992 
Higdon & Schmidt 1990 
Schmidt, Higdon, & Hueter 1988 
Matz et al. 1992 
Ogasaka et al. 1991 
Meegan et al. 1996 

(this paper) 
(this paper) 



Table 2 

Data for fitting differential peak flux distribution 



Pl 


P 2 


AN obs 


AR. (yr" 1 sr" 1 ) 


0.180 


0.202 


87 


2.45 ± 0.26 


0.202 


0.227 


83 


2.34 ± 0.26 


0.227 


0.254 


99 


2.80 ± 0.28 


0.254 


0.285 


111 


3.13 ± 0.30 


0.285 


0.320 


92 


2.60 ± 0.27 


0.320 


0.359 


96 


2.71 ± 0.28 


0.359 


0.403 


95 


2.68 ± 0.27 


0.403 


0.452 


114 


3.22 ± 0.30 


0.452 


0.507 


93 


2.62 ± 0.27 


0.507 


0.569 


79 


2.23 ± 0.25 


0.569 


0.639 


82 


2.31 ± 0.26 


0.639 


0.717 


91 


2.57 ± 0.27 


0.717 


0.804 


73 


2.06 ± 0.24 


0.804 


0.902 


83 


2.34 ± 0.26 


0.902 


1.000 


52 


1.47 ± 0.20 


1.000 


1.259 


117 


3.30 ± 0.31 


1.259 


1.584 


111 


3.13 ± 0.30 


1.584 


1.995 


104 


2.93 ± 0.29 


1.995 


2.511 


72 


2.03 ± 0.24 


2.511 


3.162 


59 


1.66 ± 0.22 


3.162 


3.981 


48 


1.35 ± 0.20 


3.981 


5.011 


40 


1.13 ± 0.18 


5.011 


6.309 


27 


0.76 ± 0.15 


6.309 


7.943 


26 


0.73 ± 0.14 


7.943 


20.00 


64 


1.80 ± 0.23 
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Table 3 

Redshifts and implied equivalent isotropic peak luminosities 



GRB 


Peak Flux 
(50-300 keV) 


z 


Implied L 
(10 51 /if 2 erg s" 1 in 30-2000 keV) 


970508 


0.97±0.05 a 


0.835 c 


0.6 ±0.1 


971214 


1.95±0.05 a 


3.418 d 


37 ± 16 


980613 


0.63±0.05 b 


1.096 6 


0.8 ±0.2 


980703 


2.39±0.06 a 


0.966 f 


2.2 ±0.4 



a Meegan et al. (1998) 

b Woods, Kippen, & Connaughton (1999) 

c Metzger et al. (1997) 

d Kulkarni et al. (1998a) 

e Djorgovski et al. (1999) 

f Djorgovski et al. (1998) 



Table 4 

Best fit parameters for monoluminous cosmological models 



p{z) 


Ro 

(/ifo Gpc -3 yr -1 sr -1 ) 


(10 51 erg s- 1 ) 


x 2 

(23 d.o.f) 


Dl (constant) 
D2 (SFR) 
D3 (AGN) 


2.0 ±0.1 
1.9 ±0.1 


40+ 06 
1.5 ±0.1 

q 1+0.4 


32.3 
64.1 
27.8 






Table 5 




Best fit parameters for cosmological models with a power-law luminosity function 


P{z) 


Ro 

(h 3 70 Gpc- 3 yr" 1 sr" 1 ) 


(10 51 h 70 2 erg s- 1 ) P 


X 2 

(22 d.o.f) 


Dl (constant) 
D2 (SFR) 
D3 (AGN) 


8.8 ± 1.3 
2.2 ±0.3 
2.1 ±0.2 


0.29±g"g| 4.6_i. 4 a 

n /1O+0.20 -, +0.3 
U - 48 -0.10 2A - 0.2 

1 44+0-53 2 fi +1 -° 


32.9 
17.9 
17.4 



a This parameter is not constrained above the the best fit value when the other parameters are free to vary. 
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Table 6 

Co-moving z = burst rate for monoluminous cosmological models 



P[Z) (h 3 70 Gpc- 3 yr- 1 ) 


GEM a 
(Myr- 1 ) 


Dl (constant) 11.9±?^ 
D2 (SFR) 0.78 ± 0.04 
D3 (AGN) 0.46 ± 0.02 


- 2 

- 0.2 

- 0.1 


a Galactic Events per Million years. 




Table 7 




Co-moving z = burst rate for cosmological models with a power-law luminosity function 


P[Z) (h 3 70 Gpc^ 3 yr- 1 ) 


GEM a 
(Myr- 1 ) 


Dl (constant) 11.1 ± 1.6 
D2 (SFR) 0.87 ±0.12 
D3 (AGN) 0.51 ± 0.05 


- 2 

- 0.2 

- 0.1 


a Galactic Events per Million years. 




Table 8 




Upper limits on fractional GRB rate due to a possible homogeneous (in Euclidean space) 


sub-population of GRBs 




p(z) 90% 


99% 


Dl (constant) 5.1% 
D2 (SFR) 6.9% 
D3 (AGN) 7.7% 


6.2% 
10.0% 
11.2% 
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